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FLUID-STRUCTURE  INTERACTION  USING  RETARDED 

POTENTIAL  AND  ABAQUS 


INTRODUCTION 

Transient  analysis  of  coupled  systems  is  in  general  a  difficult  and  demanding  problem.  Even 
with  curent  levels  of  computing,  a  good  deal  of  improvement  is  still  required  to  efficiently  analyze 
many  real  world  problems.  In  this  paper,  we  will  be  concerned  with  coupled,  fluid-structure  inter¬ 
action  problems,  in  particular,  with  submerged  structures  subjected  to  weak  shock  waves.  We  be¬ 
gin  first  with  a  brief  review  of  the  transient  analysis  of  structures  using  finite  element  methods 
(FEM).  This  is  followed  by  a  discussion  of  the  merits  of  various  approaches  for  discretizing  the 
fluid  media,  including  boundary  element  methods  (BEM)  and  the  retarded  potential  (RP)  ap¬ 
proach. 

Transient  analysis  of  linear  and  nonlinear  structures  using  the  finite  element  method  has  de¬ 
veloped  to  a  very  sophisticated  level,  see  for  instance  [1-3,22].  Finite  difference  operators  are  typ¬ 
ically  employed  for  the  time  domain  and  can  be  grouped  into  either  explicit  or  implicit  methods. 
Explicit  techniques  do  not  require  the  forming  and  factoring  of  a  global  stiffness  matrix  but  do  have 
stability  and  accuracy  concerns  which  severly  restrict  the  size  of  the  time  steps.  They  are  primarily 
applicable  to  wave  propagation  analyses,  and  have  been  very  efficiently  implemented  on  some  vec¬ 
tor  machines.  In  addition,  explicit  approaches  hold  much  promise  for  coarse  grain  parallel  process¬ 
ing  environments.  Unconditionally  stable  implicit  operators,  such  as  the  trapezoidal  and  Hilber- 
Hughes  rules  [1-3],  have  been  widely  applied  for  structural  dynamics  types  of  linear  and  nonlinear 
analyses  in  which  the  time  interval  of  concern  is  typically  much  longer  than  in  wave  propagation 
type  problems.  Implicit  methods  can  employ  much  larger  time  steps  ,  but  require  the  formuarion 
and  factoring  of  a  global  stiffness  matrix.  This  is  very  costly  and  even  prohibitive  especially  for 
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large  nonlinear  structural  problems  where  the  stiffness  matrix  must  be  repeatedly  formed  and  fac¬ 
tored.  Explicit  /implicit  methods  in  which  part  of  the  structure  is  treated  as  explicit  and  part  as  im¬ 
plicit  have  also  been  developed  but  they  too  have  similar  restrictions. 

Overall  the  transient  FEM  analysis  of  especially  large,  nonlinear  transient  structures  remains 
a  difficult  and  expensive  process.  Explicit  and  implicit  methods  both  have  advantages  and  disad¬ 
vantages  which  can  greatly  restrict  their  efficiency  -  even  on  the  largest  present  day  supercomput¬ 
ers. 


For  transient  fluid-structure  interaction  problems,  FEM  is  the  obvious  choice  for  discretization 
of  the  structure,  which  may  respond  in  a  linear  or  nonlinear  manner.  In  the  fluid  media,  which  is 
typically  approximated  as  an  acoustic  media,  several  choices  exist  for  discretization  [21].  The  han¬ 
dling  of  the  interface  between  the  fluid  and  the  structure  can  also  be  quite  involved.  Typically,  due 
to  efficency  considerations,  the  coupled  system  is  solved  in  a  staggered  fashion  [17]. 

The  treatment  of  the  fluid  media  is  motivated  by  the  fact  that  the  analyst  is  concerned  primarily 
with  the  response  of  the  structure  and  what  effect  the  fluid  has  upon  it.  Thus  additional  approxima¬ 
tions  in  the  fluid  realm  can  be  introduced  .  Two  approaches  are  evident.  The  first  is  domain  dis¬ 
cretization  with  techniques  such  as  finite  element,  finite  difference  or  finite  volume.  Of  these,  FEM 
seems  to  have  been  been  the  most  popular.  The  discretization  of  the  fluid  domain  with  finite  ele¬ 
ments,  however,  represents  a  costly  approach  that  can  also  introduce  further  approximations.  To 
achieve  a  suitable  level  of  accuracy,  a  large  number  of  elements  and  degrees  of  freedom  needs  to 
be  included  particularly  for  3D  problems,  and  typically  implicit  methods  must  be  applied.  Thus 
even  with  staggered  solutions,  this  can  be  a  very  costly  if  not  prohibitive  approach.  With  domain 
methods  such  as  FEM,  the  modeling  of  infinite  boundary  conditions  raises  significant  difficulties. 
Silent  boundary  conditions  [3]  have  been  applied  with  some  success.  However,  general  purpose 
infinite  type  finite  elements  [22]  for  transient  problems  are  not  available. 

The  second  approach  for  the  fluid  media  is  the  use  of  boundary  element  methods  (BEM)  -  see 
[4].  The  introduction  of  BEM,  while  elimimating  the  need  to  discretize  the  domain  of  the  fluid,  car¬ 
ries  with  it  an  assumption  of  homogeneity.  Thus  thermal  variations,  cavitation,  etc.  can  not  be  eas¬ 
ily  dealt  with  if  at  all.  The  retarded  potential  method  represents  an  exact  formulation  of  the  transient 
BEM  problem  [4,7,15,16,19].  However,  difficulties  especially  with  storage  have  greatly  limited  its 
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more  widespread  application.  The  doubly  asympotic  approximation  (DAA),  represents  ah  approx¬ 
imation  to  the  RP  method  and  has  been  widely  employed  for  modeling  the  transient  fluid  media  - 
see  for  instance  [5,6,9,1 1,12].  In  [  10],  DAA  was  coupled  with  a  finite  element  capability  that  mod¬ 
els  cavitation  in  the  fluid.  In  general,  however,  there  is  significant  concern  over  DAA  accuracy  es¬ 
pecially  for  highly  transient,  nonlinear  applications  of  fairly  long  durations.  DAA  correctly 
calculates  the  fluid  loading  at  early  and  late  time.  In  the  intermediate  time  region  or  for  the  radia¬ 
tion  by  the  structure  vibrations  in  the  intermediate  frequency  range,  DAA’s  accuracy  is  suspect. 
DAA2  [12]  represents  a  generalization  of  DAA  with  some  improved  accuracy  in  the  intermediate 
frequency  range. 

Overall,  the  solution  of  transient  fluid-structure  interaction  problems  is  certainly  not  at  a  fully 
satisfactory  level.  A  great  deal  of  work  remains  to  develop  efficient,  accurate  techniques  for  this 
very  important  class  of  problems. 

In  what  follows,  we  first  discuss  the  retard  potential  formulation  and  then  its  coupling  to  the 
nonlinear  finite  element  code  ABAQUS  for  the  solution  of  submerged  structures  subjected  to  weak 
shock  waves. 


FORMULATION 


Consider  a  submerged  structure  subjected  to  an  incident  pressure  wave.  We  assume  that  the 
fluid  media  can  be  approximated  as  an  ideal  compressible  medium  in  linear  wave  motion  (acoustic 
fluid).  The  incident  wave  p",cimpinges  on  the  structure  and  is  scattered  to  create  psca.  In  addition, 
the  structure’s  response  initiates  a  radiation  pressure  prad in  the  surrounding  fluid.  The  total  pres¬ 
sure  pis  then  represented  as: 


_  _  Jnc  „sca  j.  ~rad 
P  -  P  +P  +P 


(1) 


The  boundary  condition  on  the  “wetted”  surface  of  the  structure  can  be  written  as  [15] 

=  -pa  .  (2) 

where  3  is  the  unit  normal  into  the  structure;  p  is  the  fluid  density;  and  a  is  the  normal  acceleration 
of  the  structure.  The  retarded  potential  integral  is  the  solution  of  the  linear  wave  equation, 

v2P  = 

P  cW  ’ 


(3) 


subject  to  the  boundary  conditions  of  eq.  (2).  In  eq.  (3),  c  is  the  sonic  velocity  in  the  fluid. 
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A  discussion  of  the  retarded  potential  equation  and  its  discretized  approximation  as  applied  to 
submerged  structures  can  be  found  in  [15,16,19].  The  form  of  the  equation  for  the  calculation  of 
the  total  pressure  pon  the  “wetted”  surface  of  the  structure  subjected  to  a  continuous  incident  pres¬ 
sure  field  p100  is  as  follows: 

p (x,  t)  =  2 pinc (*„-(&>{ (a^jp-dS')  (4) 

+ 

s 

in  which 


t'=t- R/c 

(5) 

R  =  U-xl 

(6) 

and  where  f '  is  the  retarded  time;  x  is  the  position  of  the  observation  point  on  the  surface;  s'  is  the 
location  of  an  integration  point  on  the  surface;  R  is  the  distance  between  the  observation  point  and 
the  integration  point;  and  a'  represents  the  surface  normal  direction  at  the  integration  point  orien¬ 
tated  away  from  the  fluid. 

IMPLEMENTATION 

For  computational  purposes,  eq.  (4)  is  discretized  in  a  boundary  element  sense  and  linked  to 
the  structural  analysis  code  ABAQUS  to  form  the  program  ABAQUS-RP.  The  structure  will  be 
allowed  to  respond  in  a  linear  or  nonlinear  fashion.  The  surface  pressure  field  is  approximated  by 
subdividing  the  surface  into  K  zones  or  elements  of  constant  pressure.  These  constant  boundary 
elements  coincide  with  the  “wetted”  surface  of  the  structure  which  is  discretized  with  ABAQUS  8 
node  quadratic  shell  elements  (SR8).  Thus  the  boundary  elements  although  constant  in  the  pressure 
arc  quadratic  in  terms  of  the  geometry  [20].  Figure  1  indicates  a  8  node  ABAQUS  shell  element 
(such  as  the  S8R)  and  the  corresponding  boundary  element  divided  into  subzones.  The  pressures 
on  the  fluid  “wetted”  surface  are  discretized  in  time  to  coincide  with  the  time  steps  of  ABAQUS. 
Prior  to  conducting  the  time  history  computation,  a  matrix  of  influence  coefficents  must  be  calcu¬ 
lated  for  the  BEM  grid  which  expresses  the  relations  of  current  pressures  on  past  pressure,  and  past 
and  present  accelerations. 
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Each  zone  or  boundary  element  k  is  subdivided  into  amesh  of  subzones  /.  This  subdivision 
permits  accurate  numerical  integration  of  the  geometric  influence  coefficients  as  well  as  determi¬ 
nation  of  the  time  delayed  pressures  and  accelerations  from  each  proceeding  time  step.  For  the  pur¬ 
pose  of  measuring  the  distance  between  the  integration  point  /  and  the  observation  pointy  (thereby 
the  time  retardation),  point  j  is  located  at  the  center  of  each  boundary  element.  The  integration  point 
/  is  located  at  the  center  of  subzone  /  within  boundary  element  k.  When  the  observation  pointy  is 
within  the  zone  k  being  integrated,  R  can  go  to  zero,  and  singular  integration  techniques  must  be 
employed  -  see  [16,19]. 


The  time  step  dt  is  assumed  to  be  constant,  and  the  current  time  can  then  be  expressed  as: 


t  =  m  •  dt .  (7) 

The  retardation  time  between  points  j  and  l  (within  zone  k)  is  expressed  as  f  jkl  and  follows  from 
eq.  (S).  In  general  t'jkl  will  not  be  a  multiple  of  dt  and  may  be  expressed  as: 


*jkl  ~  njkl +  Yjkl  ’ 

where  n  is  an  integer  and  0  <  y  1 1 .  In  eq.  (4)  p,  ^  and  a  can  now  be  interpolated  as  follows: 


,m-n~  1 


dPk  .. 


(8) 


(9) 


Pk(mdt-n  =  (i -y)prn+rPk 
The  term  is  typically  expressed  by  a  three  point  backward  difference  formula  [15,16,19]  and 
we  have  employed  this  form.  However  in  [13],  Groenenboom  indicates  that  this  represents  an  ex¬ 
plicit  assumption  and  requires  for  stability  that 


dt  >  max—  .  (10) 

c 

Groenenboom  [13]  presents  an  implicit  formulation  of  the  retarded  potential  which  requires  a  ma¬ 
trix  inversion. 


The  retarded  potential  integral  eq.  (4)  can  be  discretized  and  arranged  into  the  form  [15,16,19]: 


k  _m  -  i 


(ID 


t  =  1  i  =  0 


where 


r  = 


maxRjki 

cdt 


+  1, 


/  = 


maxRikl 

- +  3  . 

cdt 


(12) 
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The  terms  and  By  depend  only  on  geometry  and  are  computed  and  stored  prior  to  calculating 
the  pressure.  Since  the  pressure  p*  on  the  right  side  of  eq.  (11)  are  not  yet  known,  a  convergence 
criteria  that  implies  that  a  disturbance  from  one  zone  cannot  affect  any  other  zone  in  less  than  one 
time  step  is  enforced,  Bqj  =  0  for  all  j*k. 

Equation  (1 1)  represents  the  response  of  the  fluid  media  on  the  “wetted”  surface  of  the  struc¬ 
ture.  Both  the  current  pressures  p*  and  current  structural  accelerations  a”  are  unknown,  and  must 
be  determined.  The  fluid  pressure  represents  a  loading  term  on  the  governing  FEM  equations  of 
the  structure.  The  coupled  equations  for  the  fluid-structure  interaction  can  be  solved  by  elimination 
of  one  field  variable,  in  a  fully  coupled  sense  (solve  both  sets  of  equations  simultaneously),  or  as 
is  done  in  this  paper  in  a  staggered  fashion  [17]  -  which  is  much  more  efficient.  With  a  staggered 
approach  the  current  pressures  p”  in  eq.  (1 1)  are  solved  using  matrix  multiplications  of  past  pres¬ 
sures  and  accelerations,  and  current  acceleration  estimates  obtained  from  the  structural  equations. 
The  solution  of  the  pressures  pj  thus  requires  storage  of  a  large  number  of  pressures  and  acceler¬ 
ations  at  previous  time  steps.  Therein  lies  the  difficulty  that  has  handicapped  the  application  of  the 
retarded  potential  approach. 

The  following  staggered  scheme  is  employed  in  this  paper: 

1 .  Apply  (p]'lc)°  at  m  =  0 

2.  Solve  the  FEM  equations  and  calculate  aj 

3.  m  =  m  +  1 

4.  From  eq.  (1 1)  calculate  pj 

5.  Apply  pj1  to  the  FEM  equations  and  calculate  aj1  * 1 

6.  If  the  maximum  time  step  is  reached,  stop. 

7.  Go  to  3. 
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We  note  that  in  using  a  staggered  approach,  a  phenomemon  known  in  control  theory  as  dead 
time  feedback  can  arise  [17].  In  [18],  Tamm  stabilized  the  solution  for  a  combined  retarded  poten¬ 
tial  -  FEM  technique  by  simply  eliminating  the  added  step  of  time  delay  in  the  structural  acceler¬ 
ation  values  relative  to  the  fluid  pressures.  This  approach  was  effective  because  of  the  rather  minor 
influence  of  die  current  accelerations  on  the  current  pressures  at  the  observation  point  j  -  see  [18] 
for  further  details.  In  this  paper,  therefore,  we  have  included  this  feature. 

APPLICATION 

A  retarded  potential  (RP)  capability  has  been  coupled  to  the  ABAQUS  program,  through  the 
DLOAO  user  written  subroutine,  to  form  ABAQUS  -  RP.  The  initial  implemention,  which  is  the 
subject  of  this  paper,  has  been  in  a  MICRO-VAX  3600  environment.  The  fluid-structure  interac¬ 
tion  problems  investigated,  therefore,  have  been  relatively  small  in  size  and  coarse  in  the  time  step¬ 
ping,  and  serve  only  to  validate  the  overall  capability.  An  important  feature  in  the  ABAQUS  -  RP 
program  is  the  Hiber-Hughes  implicit  time  operator  with  controllable  numerical  damping  [1].  Pre¬ 
vious  experience  [15,19]  has  shown  that  some  numerical  damping  has  helped  stabilize  the  fluid 
response.  In  addition,  a  variable  time  stepping  capability  has  been  included  in  the  RP  program. 
However,  it  has  proven  to  be  unstable  with  the  present  adaptive  time  stepping  algorithm  employed 
in  ABAQUS  [1].  Also  to  help  minimize  storage,  the  RP  subprogram  contains  a  capability  to  limit 
how  far  back  in  time  previous  responses  are  stored.  This  is  an  important  approximation  especially 
for  longer  time  histories  or  larger  structures,  but  our  experience  with  it  is  limited. 

To  demonstrate  the  ABAQUS-RP  program,  consider  Figure  2  in  which  an  elastic  sphere  is 
subjected  to  an  incident  step,  plane  wave  along  the  x  axis.  An  exact  solution  is  presented  in  [14], 
Figure  3  indicates  the  doubly  symmetric  ABAQUS  model  which  consists  of  54  S8R  quadratic  shell 
elements  [1]  overlaid  with  constant  pressure  boundary  elements.  Due  to  symmetry  only  1/4  of  the 
sphere  is  modeled .  The  radius  of  the  sphere  r  is  1  m,  while  the  thickness  2h  is  0.02  m.  The  prop¬ 
erties  of  the  steel  are:  Young’s  modulus  E  =  2.0684  El  1  Pa,  Poisson’s  ratio  v  =  0.3,  and  the  mass 
density  ps  =  7784.5  kg/m3 .  The  properties  of  the  surrounding  water  are:  the  speed  of  sound  c  = 
1461 .2  m/sec  and  the  mass  (tensity  pw  =  999.6  kg/m3.  The  magnitude  of  the  incident  wave  is  as¬ 
sumed  to  be  14.0E6  Pa.  The  time  step  employed  is  dr  =  0.10  r/c. 
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Figure  4  compares  the  ABAQUS  -  RP  predictions  for  the  dimensionless  displacement  history 
at  ©  =  45  and  90  degrees  with  a  closed  form  solution  [14].  With  die  above  dr,  this  model  took  ap¬ 
proximately  8.S  hours  to  go  200  time  steps  on  the  MICRO-VAX.  In  general,  the  agreement  is  only 
fair  due  probably  to  the  coarseness  of  the  time  step  and  the  grid,  and  also  to  the  relatively  high  nu¬ 
merical  damping  employed  with  the  Hilber-Hughes  operator  (a  =  0.3  [1]  was  used).  As  discussed 
earlier,  the  use  of  numerical  damping  can  be  an  important  consideration  for  the  stability  of  the  finite 
element/retarded  potential  calculations.  In  addition,  to  simulate  a  step  loading,  as  was  done  in  [15], 
a  ramp  loading  with  a  rise  time  of  t  =  .5  c/r  was  employed.  We  note  also  that  in  [15]  the  agreement 
between  NASTRAN  -  RP  results  and  the  closed  form  solution  [14]  tended  to  deteriorate  at  later 
times. 

In  an  effort  to  apply  ABAQUS-RP  to  a  nonlinear  problem,  the  sphere  is  allowed  to  yield 
(o yjeU  =  345E6  Pa  and  perfect  plasticity  is  assumed)  and  large  displacements  are  considered.  The 
magnitude  of  the  plane  wave  is  reduced  to  8.5E6  Pa  (otherwise  the  sphere  would  collapse).  The 
displacement  histories  are  are  plotted  in  Figure  5.  Comparisons  to  DAA  predictions  etc.  are  not 
available  so  this  problem  represents  only  a  demonstration  of  the  nonlinear  capability  of  the  pro¬ 
gram. 

SUMMARY  AND  FUTURE  DIRECTIONS 

In  this  paper,  we  have  presented  a  brief  discussion  concerning  transient  analysis  of  coupled 
fluid-structure  interaction  systems.  Motivated  by  the  present  state-of-the-art,  an  advanced  retarded 
potential  capability  has  been  coupled  to  the  ABAQUS  nonlinear  finite  element  program  to  produce 
ABAQUS-RP.  This  code,  which  is  executed  in  a  staggered  fashion,  is  currently  implemented  in  a 
MICRO-VAX  environment,  and  has  been  successfully  applied  to  smaller  degree  of  freedom  (dof) 
fluid-structure  interaction  problems. 

Overall,  the  RP  method  offers  significant  advantages  over  both  total  FEM  and  FEM/DAA  ap¬ 
proaches.  It  has  not  been  extensively  investigated  in  the  past  due  primarily  to  storage  requirements. 
We  have  demonstrated  that  it  can  be  fairly  efficiently  implemented  on  a  MICRO-VAX  computer 
and  applied  to  nonlinear  structures  with  a  relatively  small  number  of  dof.  Clearly,  additional  com¬ 
puter  runs  are  necessary  to  further  study  the  selection  of  time  steps,  the  use  of  numerical  damping. 
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and  the  overall  question  of  stability.  A  more  thorough  investigation  concerning  stability,  along  the 
lines  of  [17],  is  probably  in  order. 

The  remaining  challenge  is  the  implementation  of  ABAQUS-RP  in  a  supercomputer  (CRAY) 
environment.  This  will  greatly  increase  both  the  size  of  the  problems  that  could  be  efficiently  an¬ 
alyzed  and  the  length  of  the  time  history  able  to  be  considered.  In  addition,  finer  time  steps  can  be 
used.  Because  the  retarded  potential  capability  is  coupled  to  ABAQUS  through  a  user  written  sub¬ 
routine,  however,  this  will  handicap  our  ability  to  streamline  the  storage  and  retrieval  of  past  vari¬ 
ables  [7].  We  would  also  like  to  implement  a  capability  with  adaptive  time  stepping  similar  to  [8] 
and  this  may  be  difficult.  ABAQUS’ s  adaptive  time  stepping  algorithm  often  “jumps”  around  un¬ 
necessarily,  and  especially  with  expansions  can  destablize  the  fluid.  In  addition,  ABAQUS’s  im¬ 
plicit  time  stepping  approach  is  based  on  a  costly  full  Newton  formulation. 
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Structure 


S8R  Shell  Element  RP  Constant  Boundary  Element 

Figure  1. Quadratic  ABAQUS  shell  element  overlaid  by  a  constant  RP  boundary  element. 


Figure  2.  Elastic  sphere  subject  to  an  incident,  step,  plane  wave. 


X 


Figure  3.  Sphere  modeled  with  S8R  ABAQUS  shell  elements. 
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